function [ theta12Step3,out] = step3SR_theta12(namesTheta12,theta1Step1,...
    AVarTheta1Step1,theta2Step2,AVarTheta2Step2,theta12OptimalOn)

% Computing lambda and theta12Step3
numTheta12          = size(namesTheta12,1);
posTheta12Step1     = zeros(1,numTheta12);
posTheta12Step2     = zeros(1,numTheta12);
posTheta12Step2AVar = zeros(1,numTheta12);
for i=1:numTheta12
    posTheta12Step1(1,i) = find(strcmp(AVarTheta1Step1.names,namesTheta12(i)));
    posTheta12Step2(1,i) = find(strcmp(fieldnames(theta2Step2),namesTheta12(i)));
    posTheta12Step2AVar(1,i) = find(strcmp(AVarTheta2Step2.names,namesTheta12(i)));
end
theta1Step1Values = struct2array(theta1Step1)';
theta2Step2Values = struct2array(theta2Step2)';
% Note: we need this if-statement because when theta12OptimalOn = 0, then
% all covariance matrices for theta1 are NaN.
if theta12OptimalOn == 1
    lambdaRobust     = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)...
        *inv(AVarTheta1Step1.VarRobust(posTheta12Step1,posTheta12Step1) ...
        + AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar));
    
    lambdaHomo       = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)...
        *inv(AVarTheta1Step1.VarHomo(posTheta12Step1,posTheta12Step1) ...
        + AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar));
    
    lambdaHomoPooled = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)...
        *inv(AVarTheta1Step1.VarHomoPooled(posTheta12Step1,posTheta12Step1) ...
        + AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar));
    
    % Estimates of theta12Step3
    tmpRobust     = lambdaRobust*theta1Step1Values(posTheta12Step1,1) ...
        + (eye(size(lambdaRobust)) - lambdaRobust)*theta2Step2Values(posTheta12Step2,1);
    
    tmpHomo       = lambdaHomo*theta1Step1Values(posTheta12Step1,1) ...
        + (eye(size(lambdaHomo)) - lambdaHomo)*theta2Step2Values(posTheta12Step2,1);
    
    tmpHomoPooled = lambdaHomoPooled*theta1Step1Values(posTheta12Step1,1) ...
        + (eye(size(lambdaHomoPooled)) - lambdaHomoPooled)*theta2Step2Values(posTheta12Step2,1);
    
    % Computing AVar(theta12Step3)
    AVarTheta12Robust     = lambdaRobust*AVarTheta1Step1.VarRobust(posTheta12Step1,posTheta12Step1)*lambdaRobust' ...
        + (eye(size(lambdaRobust))-lambdaRobust)*AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)*(eye(size(lambdaRobust))-lambdaRobust)';
    
    AVarTheta12Homo       = lambdaHomo*AVarTheta1Step1.VarHomo(posTheta12Step1,posTheta12Step1)*lambdaHomo' ...
        + (eye(size(lambdaHomo))-lambdaHomo)*AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)*(eye(size(lambdaHomo))-lambdaHomo)';
    
    AVarTheta12HomoPooled = lambdaHomoPooled*AVarTheta1Step1.VarHomoPooled(posTheta12Step1,posTheta12Step1)*lambdaHomoPooled' ...
        + (eye(size(lambdaHomoPooled))-lambdaHomoPooled)*AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar)*(eye(size(lambdaHomoPooled))-lambdaHomoPooled)';
else
    lambdaRobust          = zeros(numTheta12,numTheta12);
    lambdaHomo            = zeros(numTheta12,numTheta12);
    lambdaHomoPooled      = zeros(numTheta12,numTheta12);
    tmpRobust             = theta2Step2Values(posTheta12Step2,1);
    tmpHomo               = theta2Step2Values(posTheta12Step2,1);
    tmpHomoPooled         = theta2Step2Values(posTheta12Step2,1);

    AVarTheta12Robust     = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar);
    AVarTheta12Homo       = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar);
    AVarTheta12HomoPooled = AVarTheta2Step2.VarRobust(posTheta12Step2AVar,posTheta12Step2AVar);
end
for i=1:numTheta12
    theta12Step3.Robust.(namesTheta12{i}) = tmpRobust(i,1);
    theta12Step3.Homo.(namesTheta12{i})   = tmpHomo(i,1);
    theta12Step3.HomoPooled.(namesTheta12{i}) = tmpHomoPooled(i,1);
end

% if no bootstrap is carried out in step 1 and/step 2, then we use the asymptotic distribution 
% to generate draws. We first fix the seeds for the random number generator
seedNum = 100;
rng(seedNum,'twister')
epsTheta1 = rand(size(AVarTheta1Step1.names,1),size(AVarTheta1Step1.names,1),AVarTheta1Step1.numBootStep1And3);
if ~isfield(AVarTheta1Step1,'RobustDraws')
    AVarTheta1Step1.RobustDraws = zeros(size(AVarTheta1Step1.names,1),AVarTheta1Step1.numBootStep1And3);
    if ~any(any(isnan(AVarTheta1Step1.VarRobust)))  
        [tmp,flag] = chol(AVarTheta1Step1.VarRobust,'lower');
        if flag ~= 0
            idx = 16;
            flagNew = 1;
            while flagNew ~= 0
                idx = idx - 1;
                [tmp,flagNew] = chol(AVarTheta1Step1.VarRobust+eye(size(AVarTheta1Step1.VarRobust,1))*10^(-idx),'lower');
            end 
        end
        for s=1:AVarTheta1Step1.numBootStep1And3
            AVarTheta1Step1.RobustDraws(:,s) = struct2array(theta1Step1)' + tmp*epsTheta1(:,s);
        end
    end
end
if ~isfield(AVarTheta1Step1,'HomoDraws')
    AVarTheta1Step1.HomoDraws = zeros(size(AVarTheta1Step1.names,1),AVarTheta1Step1.numBootStep1And3);
    if ~any(any(isnan(AVarTheta1Step1.VarHomo)))  
        [tmp,flag] = chol(AVarTheta1Step1.VarHomo,'lower');
        if flag ~= 0
            idx = 16;
            flagNew = 1;
            while flagNew ~= 0
                idx = idx - 1;
                [tmp,flagNew] = chol(AVarTheta1Step1.VarHomo+eye(size(AVarTheta1Step1.VarHomo,1))*10^(-idx),'lower');
            end 
        end
        for s=1:AVarTheta1Step1.numBootStep1And3
            AVarTheta1Step1.HomoDraws(:,s) = struct2array(theta1Step1)' + tmp*epsTheta1(:,s);
        end
    end
end
if ~isfield(AVarTheta1Step1,'HomoPooledDraws')
    AVarTheta1Step1.HomoPooledDraws = zeros(size(AVarTheta1Step1.names,1),AVarTheta1Step1.numBootStep1And3);
    if ~any(any(isnan(AVarTheta1Step1.VarHomo)))  
        [tmp,flag] = chol(AVarTheta1Step1.VarHomo,'lower');
        if flag ~= 0
            idx = 16;
            flagNew = 1;
            while flagNew ~= 0
                idx = idx - 1;
                [tmp,flagNew] = chol(AVarTheta1Step1.VarHomoPooled+eye(size(AVarTheta1Step1.VarHomoPooled,1))*10^(-idx),'lower');
            end
        end
        for s=1:AVarTheta1Step1.numBootStep1And3
            AVarTheta1Step1.HomoPooledDraws(:,s) = struct2array(theta1Step1)' + tmp*epsTheta1(:,s);
        end
    end
end
if ~isfield(AVarTheta2Step2,'Draws')
    AVarTheta2Step2.Draws = zeros(length(fieldnames(theta2Step2)),AVarTheta1Step1.numBootStep1And3);
    [tmp,flag] = chol(AVarTheta2Step2.VarRobust,'lower');
    if flag ~= 0
        idx = 16;
        flagNew = 1;
        while flagNew ~= 0
            idx = idx - 1;
            [tmp,flagNew] = chol(AVarTheta2Step2.VarRobust+eye(size(AVarTheta2Step2.VarRobust,1))*10^(-idx),'lower');
        end
    end
    epsTheta2 = rand(size(AVarTheta2Step2.names,1),size(AVarTheta2Step2.names,1),AVarTheta1Step1.numBootStep1And3);
    if length(AVarTheta2Step2.names) == length(fieldnames(theta2Step2))
        for s=1:AVarTheta1Step1.numBootStep1And3
            AVarTheta2Step2.Draws(:,s) = struct2array(theta2Step2)' + tmp*epsTheta2(:,s);
        end
    else
        % We have a reduced model in step 2
        theta2Reduced            = reduceTheta2(theta2Step2,AVarTheta2Step2.names);
        theta2Step2ValuesReduced = struct2array(theta2Reduced)';
        drawsReduced             = zeros(size(theta2Step2ValuesReduced,1),AVarTheta1Step1.numBootStep1And3);
        for s=1:AVarTheta1Step1.numBootStep1And3
            drawsReduced(:,s) = theta2Step2ValuesReduced + tmp*epsTheta2(:,s);
        end
        idx = 1;
        namesTheta2Step2 = fieldnames(theta2Step2);
        for i=1:length(namesTheta2Step2)
            if idx <= length(AVarTheta2Step2.names)
                if strcmp(namesTheta2Step2(i),AVarTheta2Step2.names(idx))
                    AVarTheta2Step2.Draws(i,:) = drawsReduced(idx,:);
                    idx = idx + 1;
                else
                    newName = namesTheta2Step2{i};
                    newName = [newName(1:end-1),'1'];
                    pos     = find(strcmp(newName,AVarTheta2Step2.names));
                    AVarTheta2Step2.Draws(i,:) = drawsReduced(pos,:);
                end
            end
        end
    end
end

% Generate the draws for teheta12 in step 3
theta12DrawsRobustStep3     = zeros(numTheta12,AVarTheta1Step1.numBootStep1And3);
theta12DrawsHomoStep3       = zeros(numTheta12,AVarTheta1Step1.numBootStep1And3);
theta12DrawsHomoPooledStep3 = zeros(numTheta12,AVarTheta1Step1.numBootStep1And3);
for s=1:AVarTheta1Step1.numBootStep1And3
    theta12DrawsRobustStep3(:,s) = lambdaRobust*AVarTheta1Step1.RobustDraws(posTheta12Step1,s) + ...
        (eye(size(lambdaRobust)) - lambdaRobust)*AVarTheta2Step2.Draws(posTheta12Step2,s); 
    theta12DrawsHomoStep3(:,s)   = lambdaHomo*AVarTheta1Step1.HomoDraws(posTheta12Step1,s) + ...
        (eye(size(lambdaHomo)) - lambdaHomo)*AVarTheta2Step2.Draws(posTheta12Step2,s); 
    theta12DrawsHomoPooledStep3(:,s) = lambdaHomoPooled*AVarTheta1Step1.HomoPooledDraws(posTheta12Step1,s) + ...
        (eye(size(lambdaHomoPooled)) - lambdaHomoPooled)*AVarTheta2Step2.Draws(posTheta12Step2,s); 
end


% The output                
out.VarRobust              = AVarTheta12Robust;
out.VarHomo                = AVarTheta12Homo;
out.VarHomoPooled          = AVarTheta12HomoPooled;
out.lambdaRobust           = lambdaRobust;
out.lambdaHomo             = lambdaHomo;
out.lambdaHomoPooled       = lambdaHomoPooled;
out.names                  = namesTheta12;
out.theta12DrawsRobust     = theta12DrawsRobustStep3;
out.theta12DrawsHomo       = theta12DrawsHomoStep3;
out.theta12DrawsHomoPooled = theta12DrawsHomoPooledStep3;
end

